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Abstract: A method is proposed for determining the masses of the new particles 
N, X, Y, Z in collider events containing a pair of effectively identical decay chains Z — > Y+ 
jet, Y — >• X + l\, X — > N + I2, where I1J2 are opposite-sign same-flavour charged leptons 
S^ and A^ is invisible. By first determining the upper edge of the dilepton invariant mass 

spectrum, we reduce the problem to a curve for each event in the 3-dimensional space of 
mass-squared differences. The region through which most curves pass then determines the 
unknown masses. A statistical approach is applied to take account of mismeasurement 
of jet and missing momenta. The method is easily visualized and rather robust against 
combinatorial ambiguities and finite detector resolution. It can be successful even for small 
event samples, since it makes full use of the kinematical information from every event. 
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1. Introduction 

One of the principal objectives of the ongoing experiments at the Large Hadron Collider 
is the discovery of new physics beyond the Standard Model. Many models of BSM physics 
predict a rich spectrum of new particles with sequential decays into chains of other new 
particles plus visible jets and leptons. Typically the endpoint of the chain is a new stable 
invisible particle that is a dark matter candidate. Important examples are the squark decay 
chain in supersymmetric models, 



q 



X2 + Q, xl^^ + F, £ t ->X?+* ± , (1-1) 



where the neutralino Xi is the lightest supersymmetric particle (LSP), and the excited 
quark decay in models with universal extra dimensions, 

q* -> Z* + q , Z* -> f ± + £ T , t± -»• 7* + ^ , (1.2) 

where the photon excitation 7* is the lightest Kaluza-Klein particle. 

If such decay chains do indeed occur at the LHC, the most urgent and challenging task 
will be to determine the masses and other properties of the new particles involved. Many 
approaches to the mass determination problem have been proposed, 1 based mainly on the 
measurement of endpoints, kinks or other features in the distributions of invariant masses 
or specially constructed observables, or on explicit solution for the unknown masses using 
multiple events. 

The present paper investigates a combination of the explicit-solution and endpoint 
methods for processes in which there are two effectively identical three-step decay chains, 



such as first- and second-generation squark pair production and decay as in (1.1). The 



method is a development of the approaches in refs. [4-6], combining kinematic fitting with 



1 For a recent review, see [1]. 

2 The explicit-solution method using pairs of events has been applied to the same class of processes in 
refs. [2,3]. 
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Figure 1: Event topology 



endpoint information to represent the possible mass solutions for each individual event as 
a curve in a three-dimensional space of mass-squared differences. For exact kinematics, 
the curves of different events all intersect at the unique correct solution point. In the 
presence of combinatorial ambiguities, measurement errors and mass variations, the region 
where the density of curves is highest gives the best estimate of the masses. Unlike pure 
endpoint or kink methods, this approach makes full use of the kinematical information 
from every event, no matter where it may lie in phase space. Determination of the four 
sparticle masses and reconstruction of the LSP momenta then appears possible even with 
small event samples. 

In Section 2 we present the general method and then in Section 3 we show results for 
a number of SUSY model points that lead to the decay chain Ql.l|). Our conclusions are 
summarized in Section 4. 



2. Method 



Consider the double decay chain in fig. HI. The 4-momenta in the upper chain should 

\2 _ „,2 

V2 



satisfy 



(pi +P2+P3+ Pif = m~ z 



(P2 +P3+P4Y 



m Y 



{P3+P4} 2 = m 2 x 

2 2 

p 4 = m N 



(2.1) 



Defining the mass-squared differences 

Mi = m% - ray > 
M 2 = my - m 2 x > 
M 3 = m\ - m% > , 

the first three equations give linear constraints on the invisible 4- momentum p^: 

- 2pi ■ p 4 = 2pi • p 2 + 2pi • p 3 + m\ - Mi = Si 



(2.2) 



-2p 2 ■ Pi = 2p 2 ■ P3 + ml - M 2 = S 2 



-2p 3 -p 4 
Similarly for the lower chain 
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M, = 5, 



(2.3) 
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We also have the missing transverse momentum constraints 
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If we make an 8-vector of the invisible 4-momenta, 

P = (pZp V 4,pI,E4,p%,p v 8 , 

then we have 



, Ea 



AP 



where A is the 8x8 matrix 
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Furthermore S may be written as 



S = BM + C 



(2.4) 



(2.5) 



(2.6) 



(2.7) 



(2.8) 



(2.9) 



where M = (Mi,M 2 ,M 3 ) is the 3- vector of mass-squared differences to be determined, 
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(2.10) 



and 



2 „x 



C = (2pi • p 2 + 2pi • p 3 + m 1 ,2p 2 ■ p 3 + m 2 , m 3 ,p ; 



miss- 



2p 5 • P6 + 2p 5 • P7 + ml 2 P6 ■ P7 + mlmlp^J . (2.11) 

Hence the solution for the invisible 4-momenta is 

p = A 1 S = DM + E (2.12) 

where D = A X B and E = A -1 C. 

The invisible 4-momenta also satisfy the quadratic constraints 

pi = Pi - Pi - Pi - Pi = m % 

pi = Pi - Pi - Pi - Pi = m% (2.13) 

where m^ is an extra unknown, independent of M. However, in the case that the terminal 
pairs of visible decay products (2,3) and (6,7) are opposite-sign same-flavour dileptons, as 



in (|1.1|) or (1.2), we can reasonably expect that the upper edge of the dilepton invariant 



mass spectrum will be measured with good accuracy. This quantity is given by 

(»D 2 = M L = {ml - m 2 x )(m 2 x - m 2 N )/m 2 x = M 2 M 3 /(M 3 + m 2 N ) (2.14) 

and hence 

m 2 N = M 3 (M 2 /M L - 1) . (2.15) 

Substituting this into eqs. (2.13), we obtain a pair of trivariate quadratic equations in 



Mi, M 2 , M 3 , whose real solutions lie on a curve in the 3-dimensional space of those variables. 
The corresponding solutions for the new particle masses are then given by eq. ( |2.15| ) and 

m 2 x = m 2 N + M 3 , m\ = m 2 x + M 2 , m 2 z = m\ + M\ . (2.16) 

The limits m 2 N > 0, m 2 z < Mjj (where My = s/4, s being the collider cm. energy squared, 
or smaller, depending on the relevant parton luminosities) imply that the solutions must 
lie within the region 

< M 3 < Mu - M L 

M L <M 2 < Mu/{\ + M 3 /M L ) 

< Mi < M v - M 2 (l + M 3 /M L ) . (2.17) 

This is a finite region with volume 

^ML^MulniMu/ML) - (3Mu - M L )(Mu - M L )] , (2.18) 

which vanishes as Ml ->0or Ml — > M\j. 
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Table 1: Mass spectrum in GeV for Snowmass point SPS la 



3. Results 



As an illustration of the method, we apply it here to the process of squark-pair production 
at the LHC (pp collisions at 14 TeV centre-of-mass energy). The SUSY mass spectrum 
and decay branching ratios are taken to be those of CMSSM point SPS la [7]. The 
corresponding masses in the decay chain (1.1) are given in Table |]. Events are generated 
using HERWIG version 6.5 [8-10]. Some of the squarks are produced directly and some 
come from gluino decay; the production mechanism affects their momentum and rapidity 
distributions but is otherwise irrelevant for our purposes. 

Third-generation squarks are excluded, as their different masses prevent a good fit 
with a single squark mass. Experimentally, this would involve vetoing events with a tagged 
6-jet. At SUSY point SPS la only left-squarks have significant branching ratios into the 
mode fll.lD and so the left-right squark mass splitting is not a problem here. The d>L — ul 
mass difference is 5.8 GeV. Therefore the assumption that the masses in the two decay 
chains are identical should be a good approximation. 

To obtain the solution curve for each event, we proceed as follows. As explained earlier, 
substitution of eq. ( 2. 15| ) into eqs. ( [2.13 ) gives a pair of trivariate quadratic equations in 
Mi, Mi-, M3. Eliminating, for example, M3 gives a quartic equation for M2 as a function of 
M\. For each real solution, M3 is given uniquely in terms M\ and M2. Thus for each M\ 
we obtain a set of 0, 2 or 4 real solution points. We divide the space of (Mi, M2, M3) into 
cells. Scanning over Mi gives a set of solution points which occupy 'hit' cells. To find all 
hit cells we also scan over M2 and M3, using permutations of the above procedure. Hits on 
already occupied cells are discarded. The resulting set of points provides an approximately 
uniform coverage of the solution curve. In general, the solution curve may consist of one 
or two closed loops or open segments with endpoints on the surface of the allowed region 



(2.17). 



Figure [2] shows the parton-level solution curves for three typical SPSla events, using 
the correct combinations of quarks and leptons in the decay chains. 3 The curves all pass 
close to the "true" mass point (TMP) 



Mi = 257040 , M 2 = 10880 , M 3 = 11233 , 



(3.1) 



all in GeV 2 , corresponding to the SUSY mass spectrum in Table |T[ The curves do not 
precisely intersect, even with exact kinematics, owing to Breit-Wigner smearing of unstable 
particle masses. However, we see that the density of solution curves is high only in the 
vicinity of the TMP (O). 



3 The two images can be merged into a three-dimensional display by directing each eye at the corre- 
sponding image. 





Figure 2: Stereoscopic views of the true parton-level solution curves for three events. The ball 
shows the true mass point. 
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Figure 3: Stereoscopic views of the parton-level solution curves for the same three events, now 
including all combinations. 



Figure |3] shows the effect of combinatorial ambiguities for the same three events, viewed 
from a different angle for clarity. Here the interchanges of near and far leptons (2 f* 3 and 
6 -<->• 7) and of quarks (1 f-> 5) are included, making eight combinations per event. Three- 
dimensional viewing reveals that incorrect combinations either have no real solutions or 
tend to give curves that do not congregate to form regions of high density. 

In the real world, the effects of parton showering, hadronization and detector reso- 
lution shift and distort the solution curves. The density of solutions around the TMP is 
reduced, and incorrect combinations may happen to produce other regions of high den- 
sity. In addition, QCD radiation produces extra jets, which increase the number of wrong 
combinations. 

To investigate these effects, we study an inclusive SUSY sample containing SUSY 
backgrounds as well as signal processes, generated using HERWIG 6.5 with initial and final 
state radiation turned on. The sample is interfaced with AcerDET 1.0 [11]. Final-state 
hadrons are formed into jets, and the momenta of jets and leptons are smeared according 
to the simulated detector resolution. 

To obtain a larger event sample with two decay chains like (|1.1| ), for this analysis 
we adopt non-CMSSM model points A, B and C, where the third-generation soft mass 
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Table 2: Parameters and mass spectra in GeV for non-CMSSM model points A, B and C. Param- 
eters common to all points are m r gen ' = 300 GeV, tan/3 = 10, sign(^) = +. 



is larger than the others, so that the branching ratio (|l.l| ) is increased by suppressing 
the X2 ~^ ^1 r=F m °de. The sparticle spectra at these points are shown in Table Q. The 
generated samples of 500,000 events correspond to about 10, 15 and 20 fb _ of integrated 
luminosity, respectively. 

The following cuts are applied in order to select signal events: 

(i) M cS = Et 1 #'* + Eti PT P,t + Et™ > 400 GeV ; 

(ii) Ef ss > max(200GeV, 0.2M eff ) ; 

(iii) At least two jets with pj, ' > 100 GeV and p ] £ ' > 50 GeV within \rj\ < 2.5 ; 

(iv) Two pairs of opposite sign same flavour leptons with px > 20 GeV and \rj\ < 3 ; 

(v) No b jet with p T > 30 GeV and \rj\ < 3. 

The b tagging efficiency is assumed to be 60%. In the cut (iv), we select not only 
opposite-flavour lepton pairs (e + e~fj, + fi~) but also the same-flavour pairs (e + e~e + e~ and 
fi + fi~fi + fi~) to have larger samples, although the latter have double the combinatorial 
background of the former. If an event contains more than two hard jets, we take the three 
hardest jets as candidates for the jets from the signal decay chains (|P|), and try all possible 
combinations. The number of combinations is 8 (16) for two candidate jets and 24 (48) for 
three with opposite (same) flavour lepton pairs. The numbers of events that survive the 
above cuts are shown in the first row in Table || together with signal/background ratios 
for each model point. The background is rather mixed, coming mainly from direct \\ 
productions associated with squarks or gauginos as well as modes containing qn — y x%3i 
b\ — > X^P an d X2 ~~ ^ Xil + l~- F° r model point C, the three-body decay x% ~^ X?^~ is 
enhanced because m^o ~ m^o + inz and turns out to be the main background. Standard 
Model background is expected to be negligible after the above selection cuts. According to 
ref. [12], the potential background comes from ti — > bbW + W~ — > 41. Based on HERWIG 
6.5 simulation of this process, we confirmed that it is indeed negligible after cuts. 

If the detector and jet properties are well understood, from the observed jet momentum, 
p ]et , we may stochastically estimate the original parton momentum, pP ar , with a gaussian 
distribution e(p par |p )et ). In this situation, we can built a confidence region in the (Mi, M2, 
M3) space [4]. For each signal event combination, i ev , a probability density function may 
be constructed as 

1 



iUM) 



N iev 



dprdpr<pr\p>r xpTw )m - om - o. (3-2) 
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326 (4.2) 
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219 (8.1) 
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Mi (True ; Best) 


231890 ; 222500 


286157 ; 282500 


316274 ; 317500 


M 2 (True ; Best) 


5624 ; 5000 


14520 ; 14200 


6815 ; 6600 


M 3 (True ; Best) 


12872 ; 11700 


10293 ; 9900 


19812 ; 18900 



Table 3: First row: number of events (signal/background) after cuts. Second row: number of 
events that contribute to the best-fit cell in the A\ 2 distribution. Third to fifth rows: true mass 
and the central value of the best- fit cell in GeV 2 . 
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Figure 4: Distribution of Ax 2 (M) for each model point at detector level. The true mass point is 
at the intersection of the three dashed lines. 



where p4, p% and m^ are the functions of M and p^ given in section ^, and Ni ev is a 
normalization factor. Given N event-combinations, log-likelihood and Ax 2 functions are 

obtained as 

N 



lnL(M) 



E 

iev 



Ww(M) 



and 



Ax 2 (M) = 2(lnL(M) r 



lnL(M)), 



(3.3) 



(3.4) 



respectively, where lnL(M) max is the maximum value of lnL(M) in the space M. 

We calculate lnL(M) approximately by the following procedure. For each event, we 
generate Monte Carlo "fake" events whose jet momenta are shifted from the original ones 
according to the probability distribution e(p par |p )et ). The parameter space M is divided 
into cells. For each cell, we count the number of fake events for which the solution curves 
go through that cell. If different combinations of the same event yield two or more curves 
passing through the same cell, we count only one. If the number of fake events is large and 
the cell size is small, this provides /i cv (M cc u) with a certain normalization. As long as we 
work with lnL(M), the normalization factor Aj cv is irrelevant, because it only shifts the 
constant term of lnL(M). We ignore cells that have /j ev (M ce u) = in our log-likelihood 
calculation, setting In /j ev (M ce u) = 0. Finally, we sum up ln/j ev (M ce u) for all combinations 
of all events. 
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Table 4: Estimated sparticle masses with their errors in GeV. 



In the following analysis, we generate 1000 Monte Carlo fake events for each event. 
For the smearing of jets and the missing transverse momentum, we use gaussian functions 
with the following standard deviations, obtained by parametrizing the AcerDET results: 



^ = -^ + 0.03, 
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for jets and 
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(3.5) 



(3.6) 



for the missing transverse momentum. We do not smear the lepton momenta because 
mismeasurement of lepton momenta is negligible compared to the jet smearing. 

Figure || shows the Ax 2 (M) distribution obtained by the above procedure for each 
model point. The cell size is AM t = 5000, AM 2 = 400, AM 3 = 600 in GeV 2 . The 
distribution has only one sharp minimum, which is close to the TMP, as can be seen in 
Table ||. Backgrounds from wrong combinations and different decay chains do not produce 
local minima at other places, and the effect of those backgrounds may be less significant 
around the true mass point. 

The second row in Table |3| shows how many different events share the best-fit cell; 
the signal/background ratios in that cell are also shown in parentheses. The ratios are 
improved significantly. For each model point the ratio is about twice that for the whole 
sample. 

In the third to fifth row of Table pL we show the central values of the best-fit cells 
compared to the TMP at each model point. As can be seen, the best-fit points are slightly 
biased towards lower masses. This may result from the following systematic errors in the 
present analysis. First, the AcerDET jets that we use are defined as massless, whereas 
the 4-momenta defined by p par = p(q) — p(x 2 ) have masses of around 10-100 GeV after 
fragmentation and hadronization. Second, we have parametrized the probability distribu- 
tions of parton momenta by gaussian functions. However, the difference between a parton 
momentum in the event record and the AcerDET jet momentum deviates slightly from a 
gaussian distribution, due to the underlying event, hadronization effects and high-p-p gluon 
emission from the original parton. A better jet algorithm with jet masses and a more 
refined parametrization will be needed to reduce these systematic errors. 

Table || shows the sparticle masses estimated from our analysis. The errors are obtained 
from la regions assuming the errors in M\, M 2 and M3 are uncorrelated, where the la 
region is defined by A% 2 < 3.53. We neglect the error from the mismeasurement of the 
dilepton endpoint because of its expected good accuracy. The la errors at Point C are 




Figure 5: One-sigma regions of 10 sub-samples, distinguished by their colours. Each sub-sample 
contains 25 events. 

accidentally small despite the small number of signal events compared to the other points. 
The sizable errors also come from the cell size, because the lcr region is almost the same 
size as the cells. More precise error estimation would require smaller cells. In addition, 
larger numbers of fake Monte Carlo events could be used to make the probability density 
smooth around the peak region. Such refinements would be justified and straightforward 
in an analysis of real data. 

The statistical approach we have adopted here is justifiable for large event samples. 
In order to see what happens in small samples, and to check the interpretation of A^ 2 , we 
divide the samples after cuts into several sub-samples, so that each of them contains 25 
events. Figure || shows the la regions of 10 sub-samples for each model point. The sub- 
samples are distinguished by their colours. As expected, the la regions are more widely 
spread compared to Figure |j. As can be seen in Figure |5]C, the gold-coloured sub-sample 
has two local minima, one of them away from the TMP. Furthermore, the la region of the 
blue-coloured sub-sample is localised away from the TMP. Those sub-samples have fewer 
signal events compared to the others. The signal/background ratio is 1.8 (2.6) for the 
blue (gold) sub-sample. We checked, however, that the la region of the blue sub-sample 
contains the TMP. Despite the small sub-sample sizes, the maximum-likelihood regions are 
still mainly localised around the TMP, and their average sizes scale as expected with the 
number of events. In our approach, unlike in endpoint and kink methods, all signal events 
contribute to the determination of the unknown masses, no matter where they may he in 
phase space, and so the method can provide meaningful information about the unknown 
masses even in rather small samples. 

4. Conclusions 

The method of mass determination presented above is simple to apply and looks promising 
for the class of processes studied here. We demonstrated the validity of the method by 
means of full simulations including detector effects. Combinatorial background, the back- 
ground from other SUSY processes and the effects of additional jets due to QCD radiation 
do not appear to be a serious problem. A statistical approach is applicable to deal with 
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jet momentum mismeasurement. We constructed an effective A% 2 variable which allows 
a rather precise determination of the unknown masses with controlled statistical errors. 
There are identified systematic errors, leading to a bias towards lower masses, which could 
be reduced with an appropriate jet algorithm and improved parametrization of jet mo- 
mentum smearing. The method can be applied successfully even to small event samples, 
because it makes full use of the kinematical information from every event. 

Note added: The procedure adopted here for constructing an approximate likelihood 
function by generating large numbers of "fake" events is quite time-consuming. An an- 
alytical procedure similar to that outlined (but not implemented) for the exact-solution 
method in Section 6 of ref. [3] may be applicable and more efficient. We thank H. C. Cheng 
for this suggestion. 
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